//首先检查是否空间折叠
//其次处理空间折叠


///这个函数必须是单形
//也就是face必须是单形
// 必须等于2 

static inline int is_edge_legacy_to_collapse(Mesh* m,template_f* f)
{
    Node*n1=m->vv_begin(m,*(f->vertices[0]));
    for(int i=1;i<f->vertices_size;i++)
    {
        if(n1==NULL){break;}
        Node* n2=m->vv_begin(m,*(f->vertices[i]));
        Node* n=node_intersection(n1,n2);
        free_node(n1);free_node(n2);
        n1=n;
    }
    int size= node_size(n1);
    free_node(n1);
    return size;
}

// M 是rows个顺序点的坐标
static double* my_deault_compute_normal(double ** M,int rows,LB_Matrix* mat)
{
    double* re=(double*)malloc(sizeof(double)*rows);
    // LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    // lb_matrix_init_double(mat,rows-1,rows-1);
    double* data=(double*)(mat->data);
    for(int k=0;k<rows;k++)
    {
        int temp_j=0;
        for(int j=0;j<rows;j++)
        {
            if(j==k){continue;}
            for(int i=0;i<rows-1;i++)
            {
                data[i*(rows-1)+temp_j]=M[i+1][j]-M[0][j];
            }
            temp_j++;
        }
        double * det=(double*)(mat->det(mat));
        re[k]=*det;
        (k%2==1)?(re[k]=-re[k]):(re[k]=re[k]);
        free(det);
    }
    //lb_matrix_free(mat);
    return re;
}

static double* my_deault_compute_normaln(template_v** cvs,int rows,LB_Matrix* mat)
{
    double* re=(double*)malloc(sizeof(double)*rows);
   
    double* data=(double*)(mat->data);
    for(int k=0;k<rows;k++)
    {
        int temp_j=0;
        for(int j=0;j<rows;j++)
        {
            if(j==k){continue;}
            for(int i=0;i<rows-1;i++)
            {
                data[i*(rows-1)+temp_j]=cvs[i+1]->point[j]-cvs[0]->point[j];
            }
            temp_j++;
        }
        double * det=(double*)(mat->det(mat));
        re[k]=*det;
        (k%2==1)?(re[k]=-re[k]):(re[k]=re[k]);
        free(det);
    }
    return re;
}


static int my_is_vertex_in_cell(template_v* cvs,int size ,template_v*v)
{
    for(int i=0;i<size;i++)
    {
        if(cvs[i]==v){return 0;}
    }
    LB_Matrix* mat =(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,size-1,size-1);
    int cols=cvs[0]->point_size;

    auto my_get_vertexp[cvs,v](int i,int j){return i==j?(v->point):(cvs[j]->point); };
    int re=1;
    if(cols==size-1)
    {
        double* data=(double*)(mat->data);
        double* inns=(double*)malloc(sizeof(double)*size);
        memset(inns,0,sizeof(double)*size);
        for(int k=0;k<size;k++)
        {
            double *p0=my_get_vertexp(k,0);
            for(int i=0;i<cols;i++)
            {
                double*  pi=my_get_vertexp(k,i+1);
                for(int j=0;j<cols;j++)
                {
                    data[i*cols+j]=pi[j]-p0[j];
                }
            }
            double* det=mat->det(mat);
            inns[k]=*det;
            free(det);
        }
        for(int i=1;i<size;i++)
        {
            if(inns[i]*inns[0]<0)
            {
                re=0;break;
            }
        }
        free(inns);
    }
    else if(cols==size)
    {
        double **M=(double**)malloc(sizeof(double*)*size);
        double* normals=(double**)malloc(sizeof(double*)*size);
        for(int k=0;k<size;k++)
        {
            for(int i=0;i<size;i++)
            {
                M[i]=my_get_vertexp(k,i);
            }
            normals[k]=my_deault_compute_normal(M,size,mat);
        }
        for(int i=0;i<size-1;i++)
        {
            for(int j=i+1;j<size;j++)
            {
                if(inner_product(normals[i],normals[j],size)<0)
                {
                    re=0;break;
                }
            }
            if(re==0){break;}
        }
        for(int i=0;i<size;i++)
        {
            free(normals[i]);
        }
        free(M);
    }
    lb_matrix_free(mat);
    return re;
}

static int my_panduan_cell_is_legal(double**M,int rows, double* normal,int cols,
    LB_Matrix* mat)
{
    int re=1;
    if(rows-1==cols)
    {
        double* data=(double*)(mat->data);
        for(int i=0;i<cols;i++)
        {
            for(int j=0;j<cols;j++)
            {
                data[i*cols+j]=M[i+1][j]-M[0][j];
            }
        }
        double* det=mat->det(mat);
        if(*det<0){re=0;}
        free(det);
        
    }
    else if(rows==cols)
    {
        if(normal==NULL){return 1;}

        double* nor=my_deault_compute_normal(M,cols,mat);
        if(inner_product(nor,normal,cols)<0)
        {
            re=0;
        }
        SAFE_FREE(nor);
    }
    return re;
}

static int my_panduan_cell_is_legaln(template_v**cvs,int rows, double* normal,int cols,
    LB_Matrix* mat)
{
    int re=1;
    if(rows-1==cols)
    {
        double* data=(double*)(mat->data);
        for(int i=0;i<cols;i++)
        {
            for(int j=0;j<cols;j++)
            {
                data[i*cols+j]=cvs[i+1]->point[j]-cvs[0]->point[j];
            }
        }
        double* det=mat->det(mat);
        if(*det<0){re=0;}
        free(det);
        
    }
    else if(rows==cols)
    {
        if(normal==NULL){return 1;}

        double* nor=my_deault_compute_normaln(cvs,cols,mat);
        if(inner_product(nor,normal,cols)<0)
        {
            re=0;
        }
        SAFE_FREE(nor);
    }

    return re;
}


static template_v*  my_is_two_halffaces_connected(template_hf* hf1,template_hf*hf2)
{
    template_v* re=NULL;
    for(int i=0;i<hf1->vertices_size;i++)
    {
        int flag=0;
        for(int j=0;j<hf1->vertices_size;j++)
        {
            if(hf1->vertices[i]==hf2->vertices[j]){flag=1;break;}
        }
        if(flag==0)
        {
            if(re==NULL){re=hf1->vertices[i];}
            else{return NULL;}
        }
    }
    return re;
}
static Int_RB_Tree* my_get_all_halffaces_from_halfface(template_hf*hf)
{
    Int_RB_Tree* tree_hf2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_hf2);

    for(int i=0;i<hf->vertices_size;i++)
    {
        for(Node* nit=hf->vertices[i]->faces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_f* f=(template_f*)(nit->value);
            template_hf*hf1=NULL;
            if(f->halffaces[0]->cell==NULL)
            {hf1=f->halffaces[0];}
            else if(f->halffaces[1]->cell==NULL)
            {hf1=f->halffaces[1];}
            if(hf1!=NULL&&hf1!=hf){tree_hf2->insert(tree_hf2,hf1->id,hf1);}
        }
    }
    return tree_hf2;
}

static void my_chuli_zhedie_cells(Mesh* m, int rows,template_v* re_nv,double* normal)
{
    LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,rows-1,rows-1);
    template_v**cvs=(template_v**)malloc(sizeof(template_v*)*rows);
    cvs[rows-1]=re_nv;
    //int mark1=0,mark0=0;
    Int_RB_Tree* tree_hf1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_hf1);

    for(auto hfit=mesh.hf_begin(&mesh);hfit!=mesh.hf_end(&mesh);hfit++)
    {
        memmove(cvs,quote(hfit)->vertices,sizeof(template_v*)*(rows-1));
        if(my_panduan_cell_is_legaln(cvs,rows,normal,re_nv->point_size,mat)==0)
        {
            tree_hf1->insert(tree_hf1,quote(hfit)->id,quote(hfit));
        }
    }
    int flag=1;
    while(flag)
    {
        flag=0;
        for(auto it =tree_hf1->begin(tree_hf1);it.it!=NULL;it++)
        {
            template_hf* hf=(template_hf*)(it.second);
            memmove(cvs,hf->vertices,sizeof(template_v*)*hf->vertices_size);
            Int_RB_Tree* tree_hf2= my_get_all_halffaces_from_halfface(hf);
            for(auto it=tree_hf2->begin(tree_hf2);it.it!=NULL;it++ )
            {
                template_hf*hf1=(template_hf*)(it.second);
                template_v* v=my_is_two_halffaces_connected(hf,hf1);
                if(v!=NULL)
                {
                    cvs[rows-1]=v;
                    if(my_panduan_cell_is_legaln(cvs,rows,normal,re_nv->point_size,mat))
                    {
                        template_c* c=  m->create_cellv(m,cvs,rows);
                        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
                        {
                            tree_hf1->erase(tree_hf1,quote(chfit)->id);
                            template_hf*hf2=Mesh_s_opposite_halfface(quote(chfit));
                            if(hf2->id==-1&&hf2->cell==NULL)
                            {
                                template_v**temp_v=(template_v**)malloc(sizeof(template_v*)*quote(chfit)->vertices_size);
                                memmove(temp_v,quote(chfit)->vertices,sizeof(template_v*)*quote(chfit)->vertices_size);
                                template_v* temp=temp_v[0];
                                temp_v[0]=temp_v[1];temp_v[1]=temp;
                                hf2=m->create_halfface(m,temp_v,quote(chfit)->vertices_size);
                                free(temp_v);
                            }
                            memmove(cvs,hf2->vertices,sizeof(template_v*)*hf2->vertices_size);
                            cvs[rows-1]=re_nv;
                            if(my_panduan_cell_is_legaln(cvs,rows,normal,re_nv->point_size,mat)==0)
                            {
                                tree_hf1->insert(tree_hf1,hf2->id,hf2);
                            }
                        }

                        flag=1;
                        break;
                    }
                }
            }
            int_rb_tree_free(tree_hf2);
            if(flag){break;}
        }
    }

    free(cvs);
    int_rb_tree_free(tree_hf1);
    lb_matrix_free(mat);
}

static Node* my_create_will_cells(Mesh* m,
    Int_RB_Tree* tree_hf,template_v* re_nv,double* normal)
{
    if(tree_hf->size<=0){return NULL;}
    Node* re=NULL;
    Mesh mesh;
    Mesh_init(&mesh);
    mesh.dimension=m->dimension;mesh.simplex=m->simplex;mesh.manifold_require=m->manifold_require;
    
    Int_RB_Tree* tree_v=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_v);
    int rows=(template_hf*)((tree_hf->begin(tree_hf)).second) +1;
    template_c* ex_c=(template_c*)malloc(sizeof(template_c));
    Cell_init_(ex_c);
    //template_hf** hfs=(template_hf**)malloc(sizeof(template_hf*)*tree_hf->size);
    int j=0;
    Node* node_v=NULL;

    for(auto it=tree_hf->begin(tree_hf);it.it!=NULL;it++)
    {
        template_hf* hf=(template_hf*)(it.second);
        template_v**temp_v=(template_v**)malloc(sizeof(template_v*)*hf->vertices_size);
        for(int i=0;i<hf->vertices_size;i++)
        {
            template_v* v=hf->vertices[i];
            template_v* nv=(template_v*)tree_v->find(tree_v,v->id);
            if(nv==NULL)
            {
                nv=mesh.create_vertexv(&mesh,v->point,v->point_size);
                node_v=node_overlying(node_v,v);
                tree_v->insert(tree_v,v->id,nv);
            }
            temp_v[i]=nv;
        }
        template_f*f= mesh.create_facev(&mesh,temp_v,hf->vertices_size);
        template_hf*hf1= mesh.create_halfface(&mesh,f,temp_v,hf->vertices_size);
        Mesh_s_opposite_halfface(hf1)->cell=ex_c;        
        //hfs[j]=hf1;
        //tree_hf1->insert(tree_hf1,hf1->id,hf1);
        free(temp_v);
        j++;
    }
   
    my_chuli_zhedie_cells( m, rows,re_nv,normal);

    template_v** temp_v=(template_v**)malloc(sizeof(template_v*)*node_size(node_v));
    int i=0;
    for(Node* nit=node_reverse(node_v);nit!=NULL;nit=(Node*)(nit->Prev))
    {
        temp_v[i]=(template_v*)(nit->value);
        i++;
    }
    for()
    {
        
    }



    free(temp_v);

    free_node(node_v);
    free_Cell(ex_c);
    int_rb_tree_free(tree_v);
    Mesh_free(mesh);
    return re;
} 

// 针对流形，且
//
template_v*  collapse_edge(Mesh* m,template_f* f ,double* p,double* normal,int cols)
{
    if(is_edge_legacy_to_collapse(m,f)>2)
    {
        return NULL;
    }

    Int_RB_Tree* tree_c=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_c);
    for(int i=0;i<f->vertices_size;i++)
    {
        for(Node* nit=f->vertices[i]->cells;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_c* c=(template_c*)(nit->value);
            tree_c->insert(tree_c,c->id,c);
        }    
    }
    //Node* bound_hfs=is_edge_legacy_to_collapse(m,f,tree_c);
    Int_RB_Tree* tree_hf=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_hf);
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
        {
            template_hf* hf=quote(chfit);
            if(tree_hf->find(tree_hf,hf->face->id)==NULL )
            {tree_hf->insert(tree_hf,hf->face->id,hf);}
            else
            {tree_hf->erase(tree_hf,hf->face->id);}
        } 
    }
    //quchu_filter_halfface(tree_hf,f);
    for(int i=0;i<f->vertices_size;i++)
    {
        for(Node* nit=f->vertices[i]->faces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_f* f=(template_f*)(nit->value);
            tree_hf->erase(tree_hf,f->id);
        }
    }
    template_v*re=m->create_vertexv(m,p,cols);



    
    
    int_rb_tree_free(tree_hf);
    int_rb_tree_free(tree_c);
    return re;
}